
############################################
############################################
# previous: 8-VAP_CDCut ####################
############################################

###### Summary ############################################################################################
#This script merges all species for each 'merging unit' together: A merging unit represents closely related species 
#with similar venom composition and similar maximum abundances, which need to be treated as one continuum of habitat 
#suitability rather than separate, overlapping taxonomic units, ie. they are not additive to each other but rather 
#'flow' from one species to the next and represent a continuous distribution that affects victims in a similar fashion.
#merging maps are created by calculating the maximum habitat suitability across all models in the merging groups and
#show groups of species that...
#are closely related AND have similar venom properties AND ...
#are likely to have one dominant species per pixel,i.e. range overlaps don't mean that they actually co-occur but rather 
#represent uncertainty of which species is present in transition zones OR for which range overlap might be true but one 
#species is likely to greatly outcompete the others in each pixel, making the main species' possible maximum abundance 
#representative of total abundance for thatgroup in that pixel.
#Once merging units are created, the second part of the script then adds their binary or continuous distributions together 
#to create hotspot maps: maps that show diversity and cumulative abundance for different 'hotspot gorups' such as all snakes, all 
#vipers,all merging groups that contain category one species, all snakes with a certain venom activity, all snakes covered by as 
#certain polyvalent antivenom, snakes that require treatment within x hours, etc. The results are plotted as .png images.
###########################################################################################################

############################################
# next: 10_VAP_ClimatChange_Summaries ######
############################################
#other potential next scripts: 
#create per administrative unit summaries of outputs to use in snakebite incidence models
#weight hotspots by actual abundance (suitability only shows abundance relative to maximum abundance 0-1), aggressiveness, etc.
#calculate likely snakebite incidence from first principals based on abundance, human exposure, etc.
############################################

#################################################################################################################
#################################################################################################################
#################################################################################################################

#Packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(stringr)
library(sf)

#Project Info:
#Project name:
projectName<-"WHOSnakes"
#Resolution in decimal degrees:
res<- 0.01
#define selection criteria of interest for further analyses:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
critThresh<-c(1)                     #c(1,1,1,10,1,75)

#Directories:
MyDir<-paste("/home/jc217070/Maxent",sep="")
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_occs")
FinalOutDir<-paste(MyDir,"/Outputs3_Final",sep="")
HotSpotDir<-paste(FinalOutDir, "/HotSpots",sep="")
IndSppOutDir<-paste(FinalOutDir, "/Final_Individual_Spp",sep="")
MergeSppOutDir<-paste(FinalOutDir, "/Final_Merge_Groups",sep="")

PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#Species Info:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
#DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin3.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax

#Spatial Data:
worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))

bg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""))
GLOBraster<-bg
values(GLOBraster)<-ifelse(values(GLOBraster)!=-Inf,0,-Inf) #create 'empty' raster of same extent/res as FAPAR
GLOBrasterB<-GLOBraster
values(GLOBrasterB)<-ifelse(values(GLOBrasterB)==0,1,values(GLOBrasterB)) #create 'empty' raster of same extent/res as hillshade

my_col1 = rev(heat.colors(n = 100, alpha=1))
my_col2<-rgb(74,122,218, max = 255, alpha = 255)
my_col3<-rgb(74,122,218, max = 255, alpha = 100)
mycols1<-c("white","grey")
mycols2<-c("royalblue","cyan","navy","skyblue","black","gold",
           "red","darkred","darkorange","orangered","deeppink","pink",
           "purple","slateblue1","lightsalmon","burlywood4","coral3",
           "deeppink4","darkorange4","cyan4","purple4")
WHOcol <- rgb(150, 180, 250, max = 255, alpha = 0) #alpha=0 means completely transparent

#####################################################################################
#####################################################################################
#####################################################################################
#####################################################################################
#####################################################################################

# A.
#################################################################################
# make binary and continuous merged (by max suitability) rasters for this group ##
# of species (e.g. all SEA spitting Cobras) for each current & future category ##
#################################################################################
########################
groupnamesA<-c(unique(DatSum$MergeShort))
groupnames<-gsub(" ","_",gsub("\\.","",gsub("\\(","",gsub("\\)","",groupnamesA))))

# groupnamesB<-c(unique(DatSum$MergeShort[DatSum$gen_sp_subtax %in% (DatSum$gen_sp_subtax[grep("NGA", DatSum$country, fixed = TRUE)]) |
#                                           DatSum$gen_sp_subtax %in% (DatSum$gen_sp_subtax[grep("IND", DatSum$country, fixed = TRUE)])]))
# groupnamesB<-gsub(" ","_",gsub("\\.","",gsub("\\(","",gsub("\\)","",groupnamesB))))

groupnamesB<-groupnames

#####

HotSpotCats<-c("Suit","Bin")

OutCats<-c("Current","2050_Median","2090_Median","2090_Q10","2090_Q90","2050_Q10","2050_Q90")
#OutCats<-c("2050_Median","2090_Median")

########################

Sys.sleep(0.1)
print(paste(Sys.time(), "making merged species group maps", sep=" "))

#####################################################################
for(groupname in groupnames){
  if(groupname %in% groupnamesB){
    
    
    
    
    
    
    # #################################################################################
    # #make temporary temp dir for rasters in this job and set from default to this dir
    # dir.create(paste(FinalOutDir,"/",groupname,"_tmpHotSpotA1",sep=""),showWarnings=FALSE)
    # rasterOptions(tmpdir=paste(FinalOutDir,"/",groupname,"_tmpHotSpotA1",sep=""))
    # #################################################################################
    # 
    k<-which(groupnames==groupname)
    myspp<-DatSum$gen_sp_subtax[DatSum$MergeShort==groupnamesA[k]]
    Sys.sleep(0.1)
    print(paste(Sys.time(), "making merged maps for",groupname, sep=" "))
    
    for(cat in OutCats){ #for each current & future category
      
      
      for(BinOrSuit in HotSpotCats){ #for binary and continuous outputs
        
        if(!file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_MergedSum.asc",sep="")) &
           !file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_MergedSum.asc.gz",sep=""))){
          
          removeTmpFiles(h=0.5) #removes all temporary files older than 30min
          
          
          
          
          
          
          ####################################      
          if(BinOrSuit=="Suit"){
            j<- 1
          } else {
            j<- 1
          }
          ####################################      
          
          SpMax<-GLOBraster
          SpSum<-GLOBraster
          
          for (sp in myspp [j:length(myspp)]){ #for each species in this subset
            
            removeTmpFiles(h=0.5) #removes all temporary files older than 30min, then read hotspot map back in to continue process
            
            if (BinOrSuit=="Bin"){
              outfile<-paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc",sep="")
            }else{
              outfile<-paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc",sep="")
            }
            
            if (file.exists(paste(outfile,".gz", sep="")) & !(file.exists(paste(outfile)))){ #safety, make sure that output .asc exists
              
              gunzip(filename=paste(outfile,".gz", sep=""), #unzip model
                     destname=paste(outfile),
                     remove=FALSE, overwrite=TRUE)
            }
            
            if (file.exists(paste(outfile, sep=""))){ #safety, make sure that output .asc exists
              
              #print out summary of whats going on (keep track of progress: start)
              Sys.sleep(0.1)
              print(paste(Sys.time(), "adding", sp, "to merge map, j=",j, sep=" "))
              
              #unzip & read thresholded sdm
              SDM<-raster(paste(outfile)) #read model output
              SDM<-merge(SDM,GLOBraster) #merge to GLOBraster to make all rasters sae extent for adding together
              
              ############################
              #calculate maximum suitability between this and previous SDMs in this group
              if(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_Merged",j-1,".asc",sep=""))){
                SpMax<-raster(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_Merged",j-1,".asc",sep="")) #read model output for species
                SpSum<-raster(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_MergedSum",j-1,".asc",sep="")) #read model output for species
              }
              
              SpMax<-max(SpMax,SDM)
              writeRaster(SpMax, filename=paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_Merged",j,".asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
              SpSum<- SpSum+SDM
              writeRaster(SpSum, filename=paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_MergedSum",j,".asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
              
              #remove previous backup file
              if(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_Merged",j-2,".asc",sep=""))){
                file.remove(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_Merged",j-2,".asc",sep=""))
                file.remove(paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_MergedSum",j-2,".asc",sep=""))
              }
              
              if (!file.exists(paste(outfile,".gz", sep="")) & (file.exists(paste(outfile)))){ #safety, make sure that output .asc exists
                gzip(filename=paste(outfile),overwrite=TRUE, remove=TRUE)
              }
              
              if (file.exists(paste(outfile,".gz", sep="")) & (file.exists(paste(outfile)))){ #safety, make sure that output .asc exists
                file.remove(filename=paste(outfile))
              }
              
            } else { #end safety loop
              
              #print out summary of whats going on (keep track of progress: start)
              Sys.sleep(0.1)
              print(paste(Sys.time(), "can't find final output for", sp, ", not added to max sutability map", sep=" "))
              
            }
            
            j<-j+1
            
          } #end species loop
          
          
          if(file.exists(paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_Merged",j-1,".asc",sep=""))){
            HotSpot<-raster(paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_Merged",j-1,".asc",sep="")) #read model output for species
            HotSpotSum<-raster(paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_MergedSum",j-1,".asc",sep="")) #read model output for species
          }
          
          if(file.exists(paste(HotSpotDir,"/",groupname,"_",cat, BinOrSuit,"_prelim_Merged",j-2,".asc",sep=""))){
            file.remove(paste(HotSpotDir,"/",groupname,"_",cat, BinOrSuit,"_prelim_Merged",j-2,".asc",sep=""))
            file.remove(paste(HotSpotDir,"/",groupname,"_",cat, BinOrSuit,"_prelim_MergedSum",j-2,".asc",sep=""))
          }
          
          maprecs<-data.frame()
          for (sp in myspp){ #for each species in this subset
            if(file.exists(paste(OccDir,"/",sp,"_allRecs.csv",sep=""))){
              #get species records used for model:
              addrecs<-read.table(paste(OccDir,"/",sp,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
              maprecs<-rbind(maprecs,addrecs)
            }}
          
          if(length(maprecs[,1])>0){
            
            Sys.sleep(0.1)
            print(paste(Sys.time(), "making extent for modelling unit ",groupname, sep=" "))
            
            rangewidthDG<-max((max(na.omit(maprecs$long))-min(na.omit(maprecs$long))),
                              (max(na.omit(maprecs$lat))-min(na.omit(maprecs$lat))))
            rangewidthM<-rangewidthDG*100000 #turn dec deg into meters
            rangewidthM<-ifelse(rangewidthM>3000000,3000000,rangewidthM) #make buffer no bigger than this
            rangewidthM<-ifelse(rangewidthM<1000000,1000000,rangewidthM) #make buffer no smaller than this
            y<-rangewidthM
            
            #make extent object describing a box around occurrence records (margin= 4 times x and y extent of occurrence records):
            xrange<- (max(na.omit(maprecs$long)))-(min(na.omit(maprecs$long)))
            yrange<- (max(na.omit(maprecs$lat)))-(min(na.omit(maprecs$lat)))
            xcenter<- (min(na.omit(maprecs$long)))+(xrange/2)
            ycenter<- (min(na.omit(maprecs$lat)))+(yrange/2)
            rangemax<-max(xrange,yrange)
            xmin<-ifelse((xcenter-(rangemax/2)- y/100000) >= -180, (xcenter-(rangemax/2)- y/100000), -180)
            xmax<-ifelse((xcenter+(rangemax/2)+ y/100000) <=  180, (xcenter+(rangemax/2)+ y/100000),  180)
            ymin<-ifelse((ycenter-(rangemax/2)- y/100000) >=  -90, (ycenter-(rangemax/2)- y/100000),  -90)
            ymax<-ifelse((ycenter+(rangemax/2)+ y/100000) <=   90, (ycenter+(rangemax/2)+ y/100000),   90)
            myext<-extent(xmin,xmax,ymin, ymax)
          } else{
            myext<-extent(-180,180,-90,90)
          }
          
          SpMax<-crop(SpMax,myext)
          SpSum<-crop(SpSum,myext)
          
          #write final version:
          writeRaster(SpMax, filename=paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_Merged.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
          writeRaster(SpSum, filename=paste(MergeSppOutDir,"/",groupname,"_",cat,BinOrSuit,"_MergedSum.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
          removeTmpFiles(h=0.5)
          
          ########################################################################
          #remove any temporary merge-files after making all hotspot maps:
          temp.file.list<-list.files(path=paste(MergeSppOutDir),pattern=paste(groupname,"_",cat,BinOrSuit,"_prelim_Merged",sep=""),full.names=TRUE)
          for (file in temp.file.list){
            file.remove(paste(file))
          }
        }}
      #######################################################################
      
      ################################################################################
      ### plot maps for this subset ##################################################
      
      Sys.sleep(0.1)
      print(paste(Sys.time(), ", plotting hotspot map", sep=" "))
      
      
      if (file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_Merged.asc.gz", sep="")) & 
          !(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_Merged.asc",sep="")))){ #safety, make sure that output .asc exists
        gunzip(filename=paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_Merged.asc.gz", sep=""), #unzip model
               destname=paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_Merged.asc",sep=""),
               remove=FALSE, overwrite=TRUE)
      }
      
      if (file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_Merged.asc.gz", sep="")) & 
          !(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_Merged.asc",sep="")))){ #safety, make sure that output .asc exists
        gunzip(filename=paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_Merged.asc.gz", sep=""), #unzip model
               destname=paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_Merged.asc",sep=""),
               remove=FALSE, overwrite=TRUE)
      }
      
      
      
      MaxBin<-raster(paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_Merged.asc",sep="")) #read model output for species
      MaxSuit<-raster(paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_Merged.asc",sep="")) #read model output for species
      
      png(filename=paste(PlotDir,"/",groupname,"_",cat,"Merged.png",sep=""), units="in", width=10, height=14, res=500) #png smaller file size, A3 for better resolution
      par(mfrow = c(2,1), mar = c(1,1,0,0), oma=c(1,1,4,1), mgp=c(1.2,0.5,0))
      
      plot(MaxSuit, col=my_col1, box=FALSE, bty="n", axes=FALSE, legend=TRUE)
      plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
      
      for (mapunit in myspp){ #for each species in this subset
        
        species_shapeA <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",paste(mapunit)))
        
        if(length(species_shapeA)>=1){
          plot(species_shapeA[2], add=TRUE, border=my_col2, col=my_col3, lwd= 1)
        }
      }
      
      plot(MaxBin, col=mycols1, box=FALSE, bty="n", axes=FALSE, legend=FALSE)
      plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
      
      for (mapunit in myspp){ #for each species in this subset
        if(file.exists(paste(OccDir,"/",mapunit,"_allRecs.csv",sep=""))){
          #get species records used for model:
          addrecs<-read.table(paste(OccDir,"/",mapunit,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
        }
        
        is<-which(DatSum$MergeShort==groupnamesA[k])
        q<-which(DatSum[is,"gen_sp_subtax"]==mapunit)
        
        species_shapeA <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",paste(mapunit)))
        
        if(length(species_shapeA)>=1){
          plot(species_shapeA[2], add=TRUE, border=mycols2[q], col=WHOcol, lwd= 1)
        }
        
        if(file.exists(paste(OccDir,"/",mapunit,".csv",sep=""))){
          urecs<-read.table(paste(OccDir,"/",mapunit,".csv",sep=""), sep=",", header=TRUE)
          if(length(urecs[,1])>=1){
            points(urecs$lat~urecs$long, col=mycols2[q], pch=1, cex=1)
            myx<-mean(urecs$long)
            myy<-min(urecs$lat)-2
            graphics::text(x=myx, y=myy, labels = gsub("_"," ",mapunit), cex = 1, col = mycols2[q])
          }
        }
      }
      
      mtext(paste(cat," maximum suitability (top) & presence/absence (bottom) for",sep=""),
            side = 3, line = 2, outer = TRUE, cex= 1.3)
      mtext(paste("'",groupnamesA[k],"'",sep=""), 
            side = 3, line = 1, outer = TRUE, cex= 1.3)
      mtext(paste("(", length(myspp), " species)",sep=""), 
            side = 3, line = 0, outer = TRUE, cex= 1.3)
      
      dev.off()
      
      
      #zip final files
      if(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_Merged.asc",sep=""))){
        gzip(filename=paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_Merged.asc",sep=""),overwrite=TRUE, remove=TRUE)
      }
      if(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_MergedSum.asc",sep=""))){
        gzip(filename=paste(MergeSppOutDir,"/",groupname,"_",cat,"Bin_MergedSum.asc",sep=""),overwrite=TRUE, remove=TRUE)
      }
      if(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_Merged.asc",sep=""))){
        gzip(filename=paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_Merged.asc",sep=""),overwrite=TRUE, remove=TRUE)
      }
      if(file.exists(paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_MergedSum.asc",sep=""))){
        gzip(filename=paste(MergeSppOutDir,"/",groupname,"_",cat,"Suit_MergedSum.asc",sep=""),overwrite=TRUE, remove=TRUE)
      }
      Sys.sleep(0.1)
      print(paste(Sys.time(), ", done plotting merged species group map", sep=" "))
      ########################################################################
      
    }
    
  }
  
  # ######################################################
  # #delete the temp dir for this script
  # unlink(paste(FinalOutDir,"/",groupname,"_tmpHotSpotA1",sep=""), recursive = TRUE) 
  # ######################################################
}

#####################################################################################
#####################################################################################
#####################################################################################
#####################################################################################
#####################################################################################
